para la predicción de Series Temporales
Índice¶
- 1 - Introducción
- 2 - Ingestión de datos y limpieza del dataset
- 3 - Comprobación de la estacionariedad de la serie (hipótesis de raíz unitaria) - R
- 4 - División en conjunto de entrenamiento y validación
- 5 - Método I: Holt-Winters con triple suavizamiento exponencial
- 6 - Método II: Predicción utilizando la librería Prophet de Facebook
- 7 - Método III: Métodos de predicción simples. Implementación en R
- 8 - Método IV: Regresión armónica dinámica (DHR) para estacionalidades múltiples - Implementado en R
- 9 - Método V: TBATS
- 10 - Conclusiones
- 11 - Futuras líneas de trabajo
- 12 - Referencias
- Anexo: Validación cruzada anidada para series temporales
0 - Nota aclaratoria¶
Toda la narrativa y explicaciones encontradas a lo largo del trabajo están hechas en español para que su correción sea más sencilla. Pero por consistencia con mis hábitos de programación, todo el código hecho en Python o R está en inglés (variables, funciones, comentarios dentro del código).
Mi intención es la de traducir todo el contenido en castellano al inglés y así poder subir el código y el resto de las explicaciones como un kernel a Kaggle como se recomendaba en la guía del TFM. Cuando tenga todo el contenido traducido al inglés lo subiré a Kaggle con el siguiente perfil: https://www.kaggle.com/avallvaq/kernels con el título - Comparison of time series forecast models. A hybrid approach in Python and R - En este apartado incluiré la nota: This analysis is part of the final project for the Master's Degree: "Master de BigData – UNED"
1 - Introducción¶
Se trada de un subconjunto de datos numéricos, entre los que existe una relación temporal.
Las series temporales necesitan su propia metodología de análisis, ya que la componente temporal añade un orden de complejidad y autocorrelación entre los datos de una misma variable.
Igual que en datos estacionarios se buscan patrones "ocultos" en los datos, ya sea entre variables del data set o variables creadas mediante proceso de data engineering. Para las series temporales se buscan la existencia de estacionalidad, análisis de tendencias o fluctuaciones periódicas (ciclos).
Existen diversos algoritmos que se pueden utilizar para analizar series temporales. Desde regresiones lineales sencillas, para hacer análisis de tendencia más sencillos, regresión mediante árboles de decisión, Random Forest o Gradient Boosting Methods o modelos ARIMA (AutoRegressive Integrated Moving Average). También existen modelos más complejos que pueden llegar a capturar las relaciones no lineales de series temporales complejas, como pueden ser las redes neuronales LSTM (Long-short Term Memory Networks), TDNN (Time Delay Neural Networks).
1.1 - Dataset de Kaggle: Hourly Energy Consumption¶
El dataset puede descargarse de Kaggle en la siguiente URL: https://www.kaggle.com/robikscube/hourly-energy-consumption#NI_hourly.csv
Se trata de una serie temporal univariante, es decir, que sólo aporta datos históricos de una sóla variable. En este caso se trata de una estimación del consumo eléctrico en Mega Watios (MW) donado por la compañia PJM Interconection LLC. Se trata de una compañía de distribución eléctrica regional en los Estados Unidos. Se encarga de distribuir electricidad entre varios estados de la zona noreste y medio este del país (algunos de ellos, Delaware, Illinois, Indiana, Kentucky, Maryland, Michigan, New Jersey, Ohio, Pennsylvania, Tennessee, Virginia, etc.)
1.2 - Análisis¶
El trabajo se ha realizado utilizando tanto Python como R. En Python se ha llevado a cabo la mayoría del EDA aprovechando la flexibilidad y potencia de la librería pandas, la mayoríad de los gráficos se renderizan en plotly, ya que la interactividad ha resultado muy útil en el contexto del análisis de series. También se ha recurrido a matplotlib y seaborn para algunas de las visualizaciones. El análisis de la serie temporal se ha llevado a cabo con la librería statsmodel y la librería Prophet de Facebook.
R es conocido por su capacidad a la hora de hacer análisis y predicción de series temporales con los objetos ts. Para el análisis de los datos se ha optado por utilizar la librería forecast, que a su vez se apoya en otros paquetes vistos durante el curso como ggplot.
El trabajo se ha desarrollado de la siguiente forma:
- Sección 2. Ingestión de datos y limpieza (eliminación de valores duplicados, creación del índice temporal, etc)
- Sección 2.1. EDA: Búsqueda de patrones estacionales, descomposición de la serie en tendecia, estacionalidad y residuos
- Sección 3. Comprobación formal de la estacionalidad de la serie. Contraste de la hipótesis de raíz unitaria (método KPPS)
- Sección 4. División simple en conjuntos de entrenamiento y validación
- Secciones 5, 6, 7, 8 y 9. Predicciones
- Sección 10. Tabla con los valores MAPE de las predicciones
1.3 - Métodos de predicción utilizados¶
En Python:
I - Holt-Winters con estacionalidad
II - Prophet
En R:
- III - Métodos sencillos en R (media, naïve, etc)
- IV - Regresión armónica dinámica con series de Fourier
- V - TBATS
1.4 - Paquetes necesarios en Python¶
import numpy as np
import pandas as pd
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()
# import pandas_profiling
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_squared_error, mean_absolute_error
import statsmodels.api as sm
import xgboost as xgb
from fbprophet import Prophet
from fbprophet.diagnostics import cross_validation
# graphic visualization packages
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import seaborn as sns
import chart_studio.plotly as py
import plotly.graph_objects as go
import plotly.io as pio
import plotly.express as px
import plotly.figure_factory as ff
%matplotlib inline
1.5 - Paquetes necesarios en R¶
# Packages needed for the R environment:
library(tidyverse)
library(ggfortify)
library(lubridate)
library(forecast)
library(tseries)
library(urca)
2 - Ingestión de datos y limpieza del dataset¶
# Read the data into a pandas DataFrame
pjme = pd.read_csv("data/PJME_hourly.csv", parse_dates=[0], encoding='UTF-8')
pjme.Datetime = pd.to_datetime(pjme.Datetime)
pjme.rename(columns={'PJME_MW':'demand_in_MW'}, inplace=True)
# Drop possible duplicate values and order the data by date and time
pjme.drop_duplicates(subset='Datetime', keep='last', inplace=True)
pjme.sort_values(by=['Datetime'], axis=0, ascending=True, inplace=True)
# Reindex after the dataframe has been sorted
pjme.reset_index(drop=True, inplace=True)
# Now, let's check if the dataframe have some missing measurments
pjme = pjme.set_index('Datetime')
print(f'Datetime freq is set to: {pjme.index.freq}, i.e. the dataset is not continous') # Datatime index frequency is set to None
# i.e. the time series is not cointinous
# We can measure the number of missing measurements
date_range = pd.date_range(start=min(pjme.index),
end=max(pjme.index),
freq='H')
print(f'The number of missing measurements in the set is: {(len(date_range)-len(pjme))}:')
# This will append the previously missing datetimes, and create null values in our target variable
pjme = pjme.reindex(date_range)
# Fill in the blanks with values that lie on a linear curve between existing data points
pjme['demand_in_MW'].interpolate(method='linear', inplace=True)
# Now the dataset is continously populated
print(f'The pjme.index.freq is now: {pjme.index.freq}, indicating that we no longer have missing instances')
# Lets copy the Datetime index again into a columns, we are going to need for further calculations
pjme['Datetime'] = pjme.index
pjme.Datetime = pd.to_datetime(pjme.Datetime)
print("\n")
# Display the first and last rows
display(pjme.head())
display(pjme.tail())
print(f"Length of the dataset: {len(pjme):,}\n\n\
Number of N/A values:\n{pd.isna(pjme).sum()}\n\
Number of missing values:\n{pd.isnull(pjme).sum()}")
2.1 - Análisis exploratorio de los datos (EDA)¶
# Let's create some seasonal features
def create_features(df, label=None):
"""
Function to create several time-range features, day of the week, day of the year, day of month,
"""
df['date'] = df.Datetime
df['hour'] = df['date'].dt.hour
df['dayofweek'] = df['date'].dt.dayofweek
df['quarter'] = df['date'].dt.quarter
df['month'] = df['date'].dt.month
df['year'] = df['date'].dt.year
df['dayofyear'] = df['date'].dt.dayofyear
df['dayofmonth'] = df['date'].dt.day
df['weekofyear'] = df['date'].dt.weekofyear
# Create the season feature
conditions = [
(df['dayofyear'] < 79),
(df['dayofyear'] >= 356),
(df['dayofyear'] >= 79) & (df['dayofyear'] < 172),
(df['dayofyear'] >= 172) & (df['dayofyear'] < 266),
(df['dayofyear'] >= 266 & (df['dayofyear'] < 356))]
choices = [0, 0, 1,2,3]
df['season'] = np.select(conditions, choices, default=0)
X = df[['date','hour','dayofweek','month', 'quarter', 'season',
'dayofyear', 'dayofmonth', 'weekofyear', 'year']]
if label:
y = df[label]
return pd.concat([X, y], axis=1)
return X
def map_months_days_to_string(df):
"""
Simple function to rename month and days of the week, from numeric values to chars.
"""
df['dayofweek_string'] = df['dayofweek'].map({0:'Monday', 1:'Tuesday', 2:'Wednesday', 3:'Thursday',
4: 'Friday', 5:'Saturday', 6:'Sunday'})
df['season_string'] = df['season'].map({0:'winter', 1:'spring', 2:'summer', 3:'fall'})
return df[['date','hour','dayofweek','dayofweek_string', 'month', 'quarter', 'season',
'season_string', 'dayofyear', 'dayofmonth', 'weekofyear', 'year', 'demand_in_MW']]
# I use the features_target dataframes for the EDA
features_target = create_features(pjme, 'demand_in_MW')
features_target = map_months_days_to_string(features_target)
# I just mapped numerical categories to strings of the actual days of the week and seasons
# for a better interpretability of the group by aggregation computations in the visual EDA section
display(features_target.head())
display(features_target.tail())
fig = go.Figure(layout=go.Layout(title=go.layout.Title(
text=f'Power Demand (MW) over time [{min(pjme.year)} - {max(pjme.year)}]')))
fig.add_trace(go.Scatter(x=pjme.date, y=pjme.demand_in_MW, name='Demand in MW', line_color='chocolate'))
fig.update_traces(line=dict(width=0.5))
fig.update_xaxes(title_text='Time')
fig.update_yaxes(title_text='Power Consumption (MW)')
fig.update_layout(showlegend=True, width=1000)
fig.show()
# Compute the correlation matrix
corr = features_target[['demand_in_MW', 'hour', 'dayofweek', 'month', 'season']].corr()
# Generate a mask for the upper triangle
mask = np.zeros_like(corr, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True
# Set up the matplotlib figure
f, ax = plt.subplots(figsize=(8, 6))
# Generate a custom diverging colormap
colormap = sns.diverging_palette(220, 10, as_cmap=True)
# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(corr, mask=mask, cmap=colormap, center=0, vmin=-1.0, vmax=1.0,
square=True, linewidths=.5, cbar_kws={"shrink": 1}, annot=True)
Podemos ver, como la hora del día tiene una correlación positiva. Por lo que en principio el consumo subirá en horas más entrado el día. A su vez también podemos ver una ligera correlación negativa del consumo eléctrico con el día de la semana. Teniendo en cuenta que hemos codificado los fines de semana como los días 5 y 6 de la semana, nos lleva a esperar una menor demanda durante el fin de semana.
2.1.1 - En busca de patrones temporales¶
# Let's aggregate the data by hour of the day during the week
features_target_aggregated = features_target.groupby(['hour', 'dayofweek_string'], as_index=False).agg({'demand_in_MW':'median'})
#plot
fig = px.line(features_target_aggregated, x='hour', y='demand_in_MW',
color='dayofweek_string',
title='Median Hourly Power Demand per Weekday')
fig.update_traces(mode='lines+markers')
fig.update_layout(xaxis_title='Hour',
yaxis_title='Energy Demand [MW]')
fig.for_each_trace(
lambda trace: trace.update(name=trace.name.replace("dayofweek_string=", "")),
)
fig.show()
Podemos observar como los fines de semana la demanda de electricidad es alrededor de un 10% menor que entre semana. A priori sería entendible plantear la hipótesis del consumo eléctrico durante la semana con el consumo que generan los puestos de trabajo.
# Let's aggregate by hour of the day during different seasons
features_target_aggregated_by_season = features_target.groupby(['hour', 'season_string'], as_index=False)\
.agg({'demand_in_MW':'median'})
# plot
color_map_season = {'winter':'mediumpurple', 'spring': 'forestgreen', 'summer': 'gold', 'fall':'chocolate'}
fig = px.line(features_target_aggregated_by_season, x='hour', y='demand_in_MW',
color='season_string', title='Median Hourly Power Demand per Season', color_discrete_map=color_map_season)
fig.update_traces(mode='lines+markers')
fig.update_layout(xaxis_title='Hour',
yaxis_title='Energy Demand [MW]')
fig.for_each_trace(
lambda trace: trace.update(name=trace.name.replace("season_string=", "")),
)
fig.show()
Aunque el coeficiente de correlación de Pearson era cercano a cero. Era esperado encontrar una estacionalidad alta del consumo eléctrico. Vemos como la demanda está muy estrechamente ligado al verano y en menor medida al invierno.
2.1.2 - Descomposición de la Serie Temporal¶
Uno de los aspectos más interesantes de las series temporales, es el poder descomponer la serie en diferentes componentes. Una componente de soporte principal de la tendencia, una componente periódica, que es la que caracteriza los ciclos estacionales.
Como en este ejemplo la componente estacional es constante, como hemos podido ver en la Figura 1.
from statsmodels.tsa.seasonal import seasonal_decompose
# NOTE, that seasonal_decompose needs to have the dataframe with a Datetime object as index
series = features_target[['demand_in_MW']]
frequency = 24*365
# decomposing the time-series, with the frequency being 24 hours per 365 days
decomposed = seasonal_decompose(series, model='additive', freq=frequency)
def plot_decompositions(decompositions, titles, line_widths):
" Plotting the different elements constituting the time-series"
for data, title, lw in zip(decompositions, titles, line_widths):
# draw a line plot of the data with plotly express
fig = px.line(data, y='demand_in_MW',
title=title, height=300)
# adjust line width
fig.update_traces(line=dict(width=lw), line_color='chocolate')
fig.update_layout(xaxis=dict(showticklabels=False, linewidth=1),
yaxis=dict(title=''),
margin=go.layout.Margin(l=40, r=40, b=0, t=40, pad=0),)
if title=='Trend':
fig.update_yaxes(range=[10000,65000])
fig.update_layout(xaxis_title="Years")
fig.show()
# calling the function
plot_decompositions(decompositions=[decomposed.trend, decomposed.seasonal, decomposed.resid],
titles=['Trend', 'Seasonality', 'Residuals'],
line_widths=[2, 0.05, 0.5])
La razón del modelo aditivo, es porque se trata de un modelo lineal, en el que la componente de tendencia (trend) tiende a ser constante, la componente estacional (seasonal) tiene una periodicidad (amplitud y frecuencia constantes).
3 - Comprobación de la estacionariedad de la serie (hipótesis de raíz unitaria) - R¶
A priori, debido a que la tendencia de la serie temporal es visualmente muy "plana", podríamos pensar que se trata de una serie estacionaria. Pero aún así, vamos a realizar una comprobación más formal de ello utilizando el método KPSS para comprobar la hipótesis de la raíz unitaria.
Hay varios métodos a disposicón para comprobar la estacionalidad de una serie temporal: -Métodos cualitativos visuales (gráfica de los datos, gráficas de autocorrelación, estadísticos móviles, etc.) -Métodos cuantitativos como el Alternative Dickey-Fuller Test (ADF) o el método KPSS
Primer punto donde introducimos código en R. Vamos a proceder a evaluar la estacionariedad de los datos, utilizando el método KPSS incluido en el paquete forecast() de R.
# First, we import the data into the R workspace
pjme = read.csv("data/PJME_hourly.csv")
pjme[1] <- lapply(pjme[1], as.character)
pjme[,1] <- parse_date_time(pjme[,1], "ymd HMS")
pjme <- distinct(pjme,Datetime, .keep_all = TRUE)
pjme <- arrange(pjme,Datetime)
pjme <- pjme[order(pjme$Datetime),]
# Train/Test split. Train: From 1st Jan 2016 to 30th Nov 2017. Test: From 1st Dec 2017 onwards
train <- pjme %>%
filter(Datetime < ymd_hms("2017-12-01 00:00:00")) %>%
filter(Datetime > ymd_hms("2015-12-31 23:00:00"))
test <- pjme %>%
filter(Datetime > ymd_hms("2017-11-30 23:00:00"))
# Create the msts() object. Which is basically a Time Series object that accepts multiple seasonality
train_ts <- msts(train$PJME_MW, start = c(2016,0), ts.frequency = 24*365.25, seasonal.periods = c(24,24*365.25))
# Since we analyzed the data and reach to the conclusion that weekly seasonality is not as strong as daily and yearly,
# I did not included it in the seasonal.periods vector, thus c(24, 24*365.25).
train_ts %>% diff() %>% ur.kpss() %>% summary()
ndiffs(train_ts) # We should take a first difference for the data
train_ts %>% diff(lag=24) %>% ndiffs() # No seasonal difference.
Pasamos a realizar las pruebas de raíz unitaria con la función KPSS del paquete urca.
# ====================
# Stationarity check
# ====================
>>> train_ts %>% ur.kpss() %>% summary() # Non-differenced data
#######################
# KPSS Unit Root Test #
#######################
Test is of type: mu with 14 lags.
Value of test-statistic is: 2.2087
Critical value for a significance level of:
10pct 5pct 2.5pct 1pct
critical values 0.347 0.463 0.574 0.739
2.21 es un valor demasiado alto del estadístico, y por lo tanto no podemos descartar la hipótesis nula de no haber una raíz unitaria en la serie temporal. Procedo a diferenciar la serie una vez más:
>>> train_ts %>% diff() %>% ur.kpss() %>% summary()
#######################
# KPSS Unit Root Test #
#######################
Test is of type: mu with 14 lags.
Value of test-statistic is: 8e-04
Critical value for a significance level of:
10pct 5pct 2.5pct 1pct
critical values 0.347 0.463 0.574 0.739
>>> ndiffs(train_ts) # We should take ONE first difference for the data
1
Tras realizar una diferenciación, los datos son estacionarios. Este es un punto muy importante del análisis, ya que el parámetro d en cualquier modelo de la familia ARIMA se decide bajo esta premisa, el número de veces que hay que diferenciar la serie hasta estacionarizarla es el valor de d.
Hacemos lo mimso con la diferencia estaciona y vemos que el resultado es nulo, por lo que no sería necesario diferenciar para hacer estacionaria la serie temporal. D=0. Es importante también tener en cuenta que el valor de D se puede ver a simple vista ya que la estacionalidad diaria (lag=24) está muy marcada y se ve que es fácilmente predecible.
>>> train_ts %>% diff(lag=24) %>% ndiffs() # No daily seasonal difference needed
0
>>> train_ts %>% diff(lag=168) %>% ndiffs() # No weekly seasonal difference needed
0
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# Let's plot the hourle autocorrelation.
fig, ax = plt.subplots(2, figsize=(12,6))
ax[0] = plot_acf(features_target.demand_in_MW.diff().dropna(), lags=168, ax=ax[0])
ax[1] = plot_pacf(features_target.demand_in_MW.diff().dropna(), lags=168, ax=ax[1])
Vemos en la gráfica de autocorrelación como la estacionalidad diaria es muy importante, ya que cada 24 horas tenemos un pico en la gráfica de autocorrelación.
Vamos a utilizar ahora la función mstl() del paquete forecast de R. Esta función nos permite descomponer la serie en distintas estacionalidades. Esto va a permitir ver además de la componente diaria, también las componentes de variación semanales y anuales.
# ============================================
# STL decomposition for multiple seasonality
# ============================================
pjme_ts %>% mstl() %>%
autoplot() + xlab("Year")
La gráfica se lee de la siguiente forma:
- Señal en crudo
- Tendencia
- Componente estacional diaria (frecuencia m=24)
- Componente estacional semanal (frecuencia m=24x7=168)
- Componente estacional anual (frecuencia m=24x365.25=8766)
Es importante fijarse en la escala vertical de las gráficas. Especialmente para ver como la estacionalidad semanal es del orden un 50% más débil que la estacionalidad diaria y anual.
4 - División en conjunto de entrenamiento y validación¶
Para el proyecto, vamos a reducir los datos a un periodo de dos años y un trimestre. Debido a la frecuencia horaria de este data-set y de la estacionalidad diaria y anual, un periodo de dos años puede contener los suficientes datos como para que nuestros modelos realicen sus estimaciones, además, reducimos el coste computacional de los modelos de predicción.
En principio y por temas de tiempo de ejecución de los distintos modelos. He obtado por una evaluación simple, con una división de los datos en entrenamiento y test del 66% 34% de la longitud de la serie. Ver Anexo para ver más información acerca de la validación cruzada de series temporales.
def ts_train_test_split(df, init_date, end_date, cutoff_date):
"""
Function to return a simple train/test split of the given a univariate time-series, within dates provided by the user.
The date should be given as a string in ISO format YYYY-MM-DD.
"""
start = init_date
end = end_date
print(f'The first date time point is: {start}')
print(f'The last date time point is: {max(df.index)}')
range_df = len(pd.date_range(start=start,end=end, freq='Y'))
print(f"The difference is: {range_df} years")
print(f'The training range has to be {range_df*0.66} years')
CUTOFF_DATE = pd.to_datetime(cutoff_date)
train = df.loc[(df.index < CUTOFF_DATE) & (df.index > init_date)].copy(deep=True)
test = df.loc[df.index >= CUTOFF_DATE].copy(deep=True)
return train, test
train, test = ts_train_test_split(features_target, '2016-01-01 00:00:00', max(features_target.index), '2017-12-01')
display(train.tail())
display(test.head())
y_train = train['demand_in_MW']
y_test = test['demand_in_MW']
display(y_test)
# Plot Train/Test
fig = go.Figure(layout=go.Layout(title=go.layout.Title(
text=f'Power Demand (MW) over time [{min(train.year)} - {max(test.year)}]')))
fig.add_trace(go.Scatter(x=train.date, y=train.demand_in_MW, name='Training', line_color='chocolate'))
fig.add_trace(go.Scatter(x=test.date, y=test.demand_in_MW, name='Test', line_color='green'))
fig.update_traces(line=dict(width=0.5))
fig.update_xaxes(title_text='Time')
fig.update_yaxes(title_text='Power Consumption (MW)')
fig.update_layout(showlegend=True, width=1000)
fig.show()
5 - Método I: Holt-Winters con triple suavizamiento exponencial¶
El algoritmo de Holt-Winters es un algoritmo basado en medias móviles muy utilizado para la predicción de series temporales. Debido a su simplicidad e interpretabilidad, vamos a utilizarlo como algoritmo base con el que después comparar los siguientes.
La mayor limitación de este algoritmo es la capacidad de tener en cuenta sólo una estacionalidad. Nos decantamos por tomar el periodo anual, seasonal_periods=24*365
import statsmodels.api as sm
# exponential smoothing only takes into consideration patterns in the target variable
# so we discard the other features
# fit & predict
holt_winter = sm.tsa.ExponentialSmoothing(y_train, seasonal='add', seasonal_periods=24*365).fit()
y_hat_holt_winter = holt_winter.forecast(len(y_test))
# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=y_test.index, y=y_test,
mode='lines',
name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=y_hat_holt_winter.index, y=y_hat_holt_winter,
mode='lines',
name='Test - Prediction'))
# adjust layout
fig.update_traces(line=dict(width=0.5))
fig.update_layout(title='Holt-Winter Forecast of Hourly Energy Demand',
xaxis_title='Date & Time',
yaxis_title='Energy Demand [MW]')
A primera vista el modelo parece obtener una buena aproximación de las series temporales. Pero vamos a evaluarlo de calculando el Mean Absolute Percentage Error (MAPE) entre el conjunto de validación y_test y la predicción y_hat_holt_winter.
def mape(y_true, y_pred):
"""
Mean Absolute Percentage Error
"""
# convert to numpy arrays
y_true, y_pred = np.array(y_true), np.array(y_pred)
# take the percentage error
pe = (y_true - y_pred) / y_true
# take the absolute values
ape = np.abs(pe)
# quantify the performance in a single number
mape = np.mean(ape)
return f'{mape*100:.2f}%'
mape_hw = mape(y_true=y_test, y_pred=y_hat_holt_winter)
print(f'Our Holt-Winters model has a mean average percentage error (MAPE) of: {mape_hw}')
A continuación se muestran las gráficas de los resultados en una ventana de tiempo menor (semanal).
# interval length
interval = 24 * 7
# intermediary variables for readability
x_true, y_true = y_test.iloc[:interval].index, y_test.iloc[:interval]
x_pred, y_pred = y_hat_holt_winter.iloc[:interval].index, y_hat_holt_winter.iloc[:interval]
# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_true, y=y_true,
mode='lines',
name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=x_pred, y=y_pred,
mode='lines',
name='Test - Prediction'))
# adjust layout
fig.update_traces(line=dict(width=0.9))
fig.update_layout(title=f'Holt-Winter Intra-Day Forecast of First {interval} Hours of Energy Demand',
xaxis_title='Date & Time',
yaxis_title='Energy Demand [MW]')
fig.show()
# quantify accuracy
print(f'MAPE for interval of the first {interval} hours: {mape(y_true, y_pred)}')
La misma ventana de tiempo pero para la última semana de predicciones.
# interval length
interval = -24 * 7
# intermediary variables for readability
x_true, y_true = y_test.iloc[interval:].index, y_test.iloc[interval:]
x_pred, y_pred = y_hat_holt_winter.iloc[interval:].index, y_hat_holt_winter.iloc[interval:]
# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_true, y=y_true,
mode='lines',
name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=x_pred, y=y_pred,
mode='lines',
name='Test - Prediction'))
# adjust layout
fig.update_traces(line=dict(width=0.9))
fig.update_layout(title=f'Holt-Winter Intra-Day Forecast of Last {abs(interval)} Hours of Energy Demand',
xaxis_title='Date & Time (yyyy/mm/dd hh:MM)',
yaxis_title='Energy Demand [MW]')
fig.show()
# quantify accuracy
print(f'MAPE for interval of the last {abs(interval)} hours: {mape(y_true, y_pred)}')
Holt-Winters MAPE: 13.11%
1st week MAPE: 5.88%
Last week MAPE: 14.24%
5.1 - El problema de la estacionalidad múltiple¶
El método de Holt-Winters con estacionalidad no está pensado para capturar estacionalidades múltiples. Al definir la frecuencia de muestreo nos vimos obligados a dar solamente un valor (24x365). A continuación vamos a utilizar modelos mas sofisticados que permiten definir estacionalidad múltiple.
6 - Método II: Predicción utilizando la librería Prophet de Facebook¶
Enlace al white paper de Prophet: https://peerj.com/preprints/3190/
Prophet se define como un método modular de regresión. Al ser modular permite muy fácilmente añadir o quitar parámetros de manera muy intuitiva, ideal para que analistas y data scientist puedan hacer uso de él y sacarle toda su potencia con una curva de aprendizaje muy suave.
La API está disponible tanto en Python como R, y es rápido, completamente automático aunque lo suficientemente flexible para que el usuario final pueda ir ajustando los modelos en función de sus necesidades o de forma iterativa.
La ventaja que tiene Prophet respecto al Método de Holt-Winters implementado en el apartado anterior, es que es capaz de lidiar con estacionalidades múltiples, eventos de un día (Black Friday, Final del Mundial...), fiestas nacionales (como Navidad, Año Nuevo, etc), periodos vacacionales, etc.
from fbprophet import Prophet
from fbprophet.diagnostics import cross_validation
# format data for prophet model using 'ds' and 'y' as per the Python API
train_prophet = y_train.reset_index().rename(columns={'index':'ds', 'demand_in_MW':'y'})
test_prophet = y_test.reset_index().rename(columns={'index':'ds', 'demand_in_MW':'y'})
display(train_prophet.tail())
display(test_prophet.head())
# Conditions for the variable seasonality of the data
def is_spring(ds):
date = pd.to_datetime(ds)
return (date.dayofyear >= 79) & (date.dayofyear < 172)
def is_summer(ds):
date = pd.to_datetime(ds)
return (date.dayofyear >= 172) & (date.dayofyear < 266)
def is_autumn(ds):
date = pd.to_datetime(ds)
return (date.dayofyear >= 266) & (date.dayofyear < 356)
def is_winter(ds):
date = pd.to_datetime(ds)
return (date.dayofyear < 79) | (date.dayofyear >= 356)
def is_weekend(ds):
date = pd.to_datetime(ds)
return date.day_name in ('Saturday', 'Sunday')
# adding to train set
train_prophet['is_spring'] = train_prophet['ds'].apply(is_spring)
train_prophet['is_summer'] = train_prophet['ds'].apply(is_summer)
train_prophet['is_autumn'] = train_prophet['ds'].apply(is_autumn)
train_prophet['is_winter'] = train_prophet['ds'].apply(is_winter)
train_prophet['is_weekend'] = train_prophet['ds'].apply(is_weekend)
train_prophet['is_weekday'] = ~train_prophet['ds'].apply(is_weekend)
# adding to test set
test_prophet['is_spring'] = test_prophet['ds'].apply(is_spring)
test_prophet['is_summer'] = test_prophet['ds'].apply(is_summer)
test_prophet['is_autumn'] = test_prophet['ds'].apply(is_autumn)
test_prophet['is_winter'] = test_prophet['ds'].apply(is_winter)
test_prophet['is_weekend'] = test_prophet['ds'].apply(is_weekend)
test_prophet['is_weekday'] = ~test_prophet['ds'].apply(is_weekend)
display(train_prophet.tail())
display(test_prophet.head())
# instantiating the Prophet class with user defined settings
prophet = Prophet(daily_seasonality=False, weekly_seasonality=False,
yearly_seasonality=False)
# custom seasonalities to account for conditional variance
# (more extreme trends in extreme seasons)
prophet.add_seasonality(name='yearly', period=365.25, fourier_order=10)
prophet.add_seasonality(name='weekly_spring',
period=7,
fourier_order=5,
condition_name='is_spring')
prophet.add_seasonality(name='weekly_summer',
period=7,
fourier_order=5,
condition_name='is_summer')
prophet.add_seasonality(name='weekly_autumn',
period=7,
fourier_order=5,
condition_name='is_autumn')
prophet.add_seasonality(name='weekly_winter',
period=7,
fourier_order=5,
condition_name='is_winter')
prophet.add_seasonality(name='daily_spring',
period=1,
fourier_order=5,
condition_name='is_spring')
prophet.add_seasonality(name='daily_summer',
period=1,
fourier_order=5,
condition_name='is_summer')
prophet.add_seasonality(name='daily_autumn',
period=1,
fourier_order=5,
condition_name='is_autumn')
prophet.add_seasonality(name='daily_winter',
period=1,
fourier_order=5,
condition_name='is_winter')
prophet.add_seasonality(name='daily_weekend',
period=1,
fourier_order=5,
condition_name='is_weekend')
prophet.add_seasonality(name='daily_weekday',
period=1,
fourier_order=5,
condition_name='is_weekday')
# fitting the model
prophet.fit(train_prophet);
# part of the dataframe on which we want to make predictions
future = test_prophet.drop(['y'], axis=1)
# predicting values
forecast = prophet.predict(future)
# see https://github.com/facebook/prophet/issues/999 for the matplotlib_converts()
pd.plotting.register_matplotlib_converters()
# plotting the seasonality components found
prophet.plot_components(forecast)
# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=test_prophet.ds, y=test_prophet.y,
mode='lines',
name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=forecast.ds, y=forecast.yhat,
mode='lines',
name='Test - Prediction'))
# adjust layout
fig.update_traces(line=dict(width=0.5))
fig.update_layout(title='Prophet Forecast of Hourly Energy Demand',
xaxis_title='Date & Time (yyyy/mm/dd hh:MM)',
yaxis_title='Energy Demand [MW]')
fig.show()
# quantify accuracy
print(f'MAPE for Prophet\'s predictions: {mape(test_prophet.y, forecast.yhat)}')
# interval length
interval = 24 * 7
# intermediary variables for readability
x_true, y_true = test_prophet.iloc[:interval].ds, test_prophet.iloc[:interval].y
x_pred, y_pred = forecast.iloc[:interval].ds, forecast.iloc[:interval].yhat
# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_true, y=y_true,
mode='lines',
name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=x_pred, y=y_pred,
mode='lines',
name='Test - Prediction'))
# adjust layout
fig.update_traces(line=dict(width=0.9))
fig.update_layout(title=f'Prophet Intra-Day Forecast of First {interval} Hours of Energy Demand',
xaxis_title='Date & Time (yyyy/mm/dd hh:MM)',
yaxis_title='Energy Demand [MW]')
fig.show()
# quantify accuracy
print(f'MAPE for interval of the first {interval} hours: {mape(y_true, y_pred)}')
# interval length
interval = -24 * 7
# intermediary variables for readability
x_true, y_true = test_prophet.iloc[interval:].ds, test_prophet.iloc[interval:].y
x_pred, y_pred = forecast.iloc[interval:].ds, forecast.iloc[interval:].yhat
# create figure
fig = go.Figure()
fig.add_trace(go.Scatter(x=x_true, y=y_true,
mode='lines',
name='Test - Ground Truth'))
fig.add_trace(go.Scatter(x=x_pred, y=y_pred,
mode='lines',
name='Test - Prediction'))
# adjust layout
fig.update_traces(line=dict(width=0.9))
fig.update_layout(title=f'Prophet Intra-Day Forecast of Last {abs(interval)} Hours of Energy Demand',
xaxis_title='Date & Time (yyyy/mm/dd hh:MM)',
yaxis_title='Energy Demand [MW]')
fig.show()
# quantify accuracy
print(f'MAPE for interval of the last {abs(interval)} hours: {mape(y_true, y_pred)}')
Prophet's MAPE: 9.06%
1st week forecast MAPE: 6.71%
Last wekk forecast MAPE: 11.89%
7 - Método III: Métodos de predicción simples. Implementación en R¶
Para la parte del trabajo realizada en R me he decantado por utilizar el paquete forecast de Rob Hyndman https://github.com/robjhyndman/forecast. La documentación está muy detallada, y además es posible encontrar numerosos ejemplos en su libro Forecasting: Principles and Practice.
Como ya vimos en la sección 3 de de las pruebas de raíz unitaria. El objeto base del paquete forecast es la clase ts incluida en el paquete básico de R. Forecast generaliza ts introduciendo la clase msts. El primero es un objeto time-series con frecuencia regular y estacionalidad única. Mientras que msts aunque también es un objeto de frecuencia regular (no permite fechas ni valores vacíos como un objeto xts), permite pasarle un vector con las frecuencias más importantes (diarias, semanales, anuales, etc). Esta flexibilidad veremos que puede marcar la diferencia a la hora de realizar predicciones en función de la serie temporal a tratar.
Para las siguientes secciones del trabajo he seguido un proceso incremental. Primero he analizado la serie con los métodos más sencillos (aunque bastante robustos y eficientes) y he ido aumentando la complejidad de los modelos utilizados (Regresión armónica dinámica con series de Fourier y finalmente TBATS).
En los siguientes puntos de la sección 7, mostraré los resultados obtenidos con estos métodos más sencillos.
7.1 - Mean method¶
Este método consiste simplemente en dar a todos los valores futuros, el valor medio de la serie.
>>> pjme_mean <- meanf(train_ts, h=24*245)
# Plot
>>> autoplot(train_ts) +
autolayer(test_ts, series="Test") +
autolayer(pjme_mean, series="Mean", PI=FALSE) +
ggtitle("Power Demand (MW) over time [2016 - 2018]") +
xlab("Year") + ylab("Demand in MW") +
guides(colour=guide_legend(title="Forecast"))
# MAPE (%) computation
>>> acc1 <- round(accuracy(pjme_mean, test_ts),2)[2,5]
>>> sprintf("Mean method MAPE: %s %%", acc1)
"Mean method MAPE: 14.7 %"
7.2 - Naïve method¶
En este caso, todas las predicciones reciben el valor del último dato disponible.
>>> pjme_naive <- naive(train_ts, h=24*245)
# Plot
>>> autoplot(train_ts) +
autolayer(test_ts, series="Test") +
autolayer(pjme_naive, series="Naïve", PI=FALSE) +
ggtitle("Power Demand (MW) over time [2016 - 2018]") +
xlab("Year") + ylab("Demand in MW") +
guides(colour=guide_legend(title="Forecast"))
# MAPE (%) computation
>>> acc2 <- round(accuracy(pjme_naive, test_ts),2)[2,5]
>>> sprintf("Naive method MAPE: %s %%", acc2)
"Naive method MAPE: 15.18 %"
7.3 - Naïve with drift method¶
En este método se permite aumentar o disminuir el valor de las predicciones en el tiempo. Es como si trazasemos una línea entre el primer punto histórico de la serie y el último y extrapolásemos con la misma pendiente a valores futuros.
>>> pjme_drift <- rwf(train_ts, h=24*245, drift=TRUE)
# Plot
>>> autoplot(train_ts) +
autolayer(test_ts, series="Test") +
autolayer(pjme_drift, series="Drift", PI=FALSE) +
ggtitle("Power Demand (MW) over time [2016 - 2018]") +
xlab("Year") + ylab("Demand in MW") +
guides(colour=guide_legend(title="Forecast"))
# MAPE (%) computation
>>> acc3 <- round(accuracy(pjme_drift, test_ts),2)[2,5]
>>> sprintf("Naive method with drift MAPE: %s %%", acc3)
"Naive method with drift MAPE: 15.04 %"
7.4 - Seasonal Naïve method¶
Los visualización de las gráficas dadas por este método parecen más complicadas de lo que en realidad son. Pero simplemente los puntos de las predicciones son los valores en estacionalidades anteriores. Es decir, si tenemos un valor y_t+1 = y_t en la estacionalidad pasada. Por ejemplo los puntos en Diciembre de 2019 tendrán el valor de puntos en Diciembre de 2018.
>>> pjme_s_naive <- snaive(train_ts, h=24*245)
# Plot
>>> autoplot(train_ts) +
autolayer(test_ts, series="Test") +
autolayer(pjme_s_naive, series="Seasonal naïve", PI=FALSE, alpha=0.7) +
ggtitle("Power Demand (MW) over time [2016 - 2018]") +
xlab("Year") + ylab("Demand in MW") +
guides(colour=guide_legend(title="Forecast"))
# MAPE (%) computation
>>> acc4 <- round(accuracy(pjme_s_naive, test_ts),2)[2,5]
>>> sprintf("Seasonal Naive method MAPE: %s %%", acc4)
"Seasonal naïve method MAPE: 13.12 %"
8 - Método IV: Regresión armónica dinámica (DHR) para estacionalidades múltiples - Implementado en R¶
Cuando hay estacionalidades muy largas, una regresión armónica dinámica utilizando series de Fourier para modelizar los regresores externos es más indicado que métodos como ARIMA o ETS. En el dataset escogido para este trabajo, al tener datos por hora, estamos hablando de que cualquier estacionalidad presente va a ser muy alta, m=24 para estacionalidad diaria, m=168 para una estacionalidad semanal, m=8766 para estacionalidades anuales. La principal diferencia es que los métodos de la familia ARIMA o ETS están concebidos para estacionalidades más bajas (mensuales, m=12 o trimestrales m=4). La principal razón es de coste computacional, ya que estos métodos necesitan calcular m-1 parámetros, con el consiguiente coste computacional y de memoria.
En los métodos de regresión armónica, la componente estacional se calcula añadiendo transformadas de Fourier en las distintas frecuencias de cada estacionalidad modelada. El movimiento periódico de corta duración de la serie se modeliza con un modelo ARMA (modelo autoregresivo de media móvil).
Las principales ventajas que ofrece este método son:
- Estacionalidad de cualquier duración (para diferentes estacionalidad, se añaden distintos términos de Fourier con distintas frecuencias).
- La suavidad del patrón estacional puede controlarse con K, el número de pares de cosenos y senos en cada Transformada de Fourier para cada frecuencia. Un valor más bajo de K implica un modelado más suave.
- El movimiento periódico a corto plazo se modela con un simple modelo ARMA.
8.1 - Implementación en R¶
Para implementar este método, me he decantado por utilizar la librería forecast() de Rob Hyndman. La forma de implementarlo es la siguiente:
Primero utilizamos la función auto.arima().
fit <- auto.arima(train_ts, seasonal=TRUE, lambda = 0, xreg = fourier(train_ts, K=c(12,2)))
xreg, es el parámetro más importante porque es el que se encarga de modelizar las estacionalidades y añadirlas al modelo ARMA.
La componente ARMA del modelo lo selecciona de forma automática utilizando como criterio la minimización del AICc (bias corrected Akaike information criterion).
>>> fit # Output resumen del modelo:
Series: train_ts
Regression with ARIMA(5,1,5) errors
Box Cox transformation: lambda= 0
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3 ma4 ma5
1.6553 -1.2152 1.4254 -1.2693 0.3028 -0.6828 -0.0202 -0.8966 0.289 0.4407
s.e. 0.0142 0.0229 0.0197 0.0202 0.0126 0.0134 0.0142 0.0070 0.013 0.0096
drift S1-24 C1-24 S2-24 C2-24 S3-24 C3-24 S4-24 C4-24 S5-24
0e+00 -0.1268 -0.0817 -0.0644 -0.0042 -0.0041 0.0026 0.0044 -0.0036 -0.0012
s.e. 1e-04 0.0042 0.0043 0.0007 0.0007 0.0004 0.0004 0.0003 0.0003 0.0002
C5-24 S6-24 C6-24 S7-24 C7-24 S8-24 C8-24 S9-24 C9-24 S10-24 C10-24
-0.0043 -1e-03 -7e-04 7e-04 7e-04 0 1e-04 -1e-04 0 -1e-04 0
s.e. 0.0002 1e-04 1e-04 1e-04 1e-04 0 0e+00 0e+00 0 0e+00 0
S11-24 C11-24 C12-24 S1-168 C1-168 S2-168 C2-168 S1-8766 C1-8766 S2-8766
0 0 0 -0.0439 0.0174 0.0211 0.0233 0.0191 -0.0330 0.1163
s.e. 0 0 0 0.0053 0.0053 0.0025 0.0026 0.3038 0.2747 0.1426
C2-8766 S3-8766 C3-8766 S4-8766 C4-8766
0.0976 -0.0064 -0.0053 0.0058 0.0081
s.e. 0.1387 0.0936 0.0935 0.0695 0.0706
sigma^2 estimated as 0.0002021: log likelihood=47630.77
AIC=-95167.55 AICc=-95167.28 BIC=-94804.29
Una vez tenemos el modelo, el siguiente paso es hacer una gráfica de las predicciones frente a los valores de test, y comprobar que los residuos tengan una apariencia de ruido blanco:
>>> fit <- auto.arima(train_ts, seasonal=TRUE, lambda = 0, xreg = fourier(train_ts, K=c(12,2,4)))
# Plot:
fit %>%
forecast(xreg=fourier(train_ts, K=c(12,2,4), h=24*245), level = 99) %>%
autoplot(include=945*24, PI=FALSE) +
autolayer(test_ts, series="Test", alpha=0.8) +
ylab("Demand in MW") + xlab("Date")
checkresiduals(fit)
"Dynamic Regression with Fourier terms MAPE: 1.04 %"
A la hora de interpretar los residuos de este modelo, como con otras técnicas de regresión, lo principal es que los residuos tengan una distribución, normal, aleatoria o de ruido blanco como se conoce comúnmente en el campo de procesado de señales.
9 - Método V: TBATS¶
Alternativa a la Regresión armónica dinámica con términos de Fourier. Está desarrollada por De Livera, Hyndman, & Snyder (https://robjhyndman.com/publications/complex-seasonality/).
La función tbats() del paquete forecast de R está completamente automatizada, aunque el usuario tiene cierta libertad para seleccionar algunos parámetros de entrada. La elección de estos parámetros reduce el número de combinaciones posible y por lo tanto disminuye el número total de modelos a elegir, por lo que la búsqueda del mejor modelo será más eficiente.
Cabe destacar la posibilidad de pasar a la función el parámetro use.parallel=TRUE. Con esta opción seleccionada tbats() aprovechará la capacidad multi-hilo del procesador, con la consecuente reducción del tiempo de cálculo.
>>> train_ts %>%
tbats(use.trend = FALSE, use.arma.errors = TRUE,
use.parallel = TRUE) -> fit2
>>> fit2
TBATS(1, {5,4}, -, {<24,1>, <168,1>, <8766,1>})
Call: tbats(y = ., use.trend = FALSE, use.arma.errors = TRUE, use.parallel = TRUE)
Parameters
Alpha: 0.1098577
Gamma-1 Values: 0.000275106 -4.054186e-05 -0.0003767173
Gamma-2 Values: 0.0003092811 -0.0002484964 -0.0001945771
AR coefficients: 2.419955 -2.72113 1.520097 -0.248297 -0.125298
MA coefficients: -0.445064 0.253288 0.504738 0.15711
# Forecast and Plot:
>>> fc2 <- forecast(fit2, h=24*245, level = 80)
>>> autoplot(train_ts) +
autolayer(test_ts, series="Test") +
autolayer(fc2, series="TBATS", PI=FALSE, alpha=0.7) +
ggtitle("Power Demand (MW) over time [2016 - 2018]") +
xlab("Year") + ylab("Demand in MW") +
guides(colour=guide_legend(title="Forecast"))
"TBATS MAPE: 1.22 %"
Vemos como el modelo TBATS ha tenido peor rendimiento que el modelo DHR. Como se puede ver en las imágenes no ha conseguido capturar de forma tan exacta la variación estacional de mayor amplitud.
9.1 - Resumen de los resultados¶
| Model | MAPE (%) | Language |
|---|---|---|
| Holt-Winters | 13.11 | Python |
| Prophet | 9.06 | Python |
| Mean | 14.7 | R |
| Naïve | 15.18 | R |
| Naïve with drift | 15.04 | R |
| Seasonal naïve | 13.12 | R |
| Dynamic harmonic Reg. + Fourier | 1.04 | R |
| TBATS | 1.22 | R |
10 - Conclusiones¶
Las conclusiones podría resumirlas en los siguientes puntos:
- El modelo de Holt-Winters, es un modelo que no está pensado para tratar series temporales con estacionalidad múltiple, aún así ha parecido sacar una solución bastante robusta.
- La librería Prophet de Facebook, ofrece de manera sencilla una doble API (en Python y R) y muy potente un modelo de predicción con variables de estacionalidad por Transformadas de Fourier. El modelo ha resultado ser muy robusto, sin apenas haber experimentado con los parámetros, y además ofrece una gran flexibilidad, añadir un gran número de estacionalidades de diferentes frecuencias, funciona con datos intra día automáticamente si se le inyecta la serie temporal en forma de dataframe pandas con un índice de tipo Datetime.
- El paquete forecast de R otorga al usuario la posibilidad de trabajar de forma muy potente con gran cantidad de modelos, transformacioens y utilidades para analizar y predecir series temporales. Al igual que con la librería Prophet, la experiencia y el tiempo son factores clave para seleccionar los parámetros adecuados para cada serie temporal.
- Se ha podido observar durante la realización de este trabajo la gran cantidad de factores y parámetros que hay que tener en cuenta para ajustar los modelos de predicción a una serie temporal. Por esta razón modelos como el de Facebook, cuyos valores han sido elegidos en base a ejemplos en la documentación, son con un gran alto de probabilidad subóptimos. Es muy notable la diferencia con los modelos automáticos implementados desde el paquete forecast de R. La selección automática a través de la minimización del AICc ha dado lugar a la obtención de modelos más indicados para la serie temporal que se estaba analizando.
11 - Futuras líneas de trabajo¶
- En este proyecto se ha considerado una serie temporal univariante, sin ningún tipo de serie adicional variable complementaria ni covariantes (por ejemplo variables climatológicas como la temperatura, horas de sol, etc). Este es uno de los puntos débiles de este análisis, porque es precisamente la fuerte correlación entre el tiempo y el consumo eléctrico uno de los puntos fuertes de los modelos de demanda de electricidad. Este comportamiento se ha visto reflejado durante el EDA, cuando los patrones de gasto eléctrico diferían mucho entre las estaciones del año.
- Los fines de semana, los periodos estivales (cuando es mayor la cantidad de gente de vacaciones), los eventos de un día afectan en gran medida al consumo. Así como la propia localización de los datos. No es lo mismo tomar datos de demanda eléctrica en una urbe, que en una zona rural, o en una zona turística donde se esperaran aún mayor varianza entre los periodos de vacaciones y el resto del año. Por todos estos factores es muy común la presencia de outliers presentes en los datos en forma de picos de gasto o días de muy poco consumo. Este es sin duda uno de los puntos fuerte de la Prohet de Facebook, que es la facilidad para modelizar este tipo de casos especiales.
- Automatizar la selección del hiperparámetro K para los modelos en los que se modelizan regresores externos de estacionalidad mediante series de Fourier. Implementando un pequeño Grid Search entre diferentes valores de K, y tomando como condición una minimización del valor del AICc podría dar lugar a mejores predicciones.
- Por último, pienso que sería interesante automatizar una validación cruzada con diferentes rangos de la serie temporal. De esta forma se podrían ejectuar diferentes modelos con diferentes rangos de datos y obtener finalmente un modelo más robusto.
La implementación de la validación cruzada podría partir de este ciclo junto con el fragmento de código que hay en el Anexo:
for train_index, test_index in tscv.split(X):
print("TRAIN:", train_index, "TEST:", test_index)
X_train, X_test = X.iloc[train_index], X.iloc[test_index] # pandas doesn't support directly integer based slicing, so that .iloc is mandatory
y_train, y_test = y.iloc[train_index], y.iloc[test_index]
X_train, y_train = model_forecast.fit()
y_pred = model_forecast.predict(X_test) # in literature is usually seen y_pred as y_hat
cross_validaiton_result(y_test, y_pred)
12 - Referencias¶
- Hyndman, R.J., & Athanasopoulos, G. (2018) Forecasting: principles and practice, 2nd edition, OTexts: Melbourne, Australia. https://otexts.com/fpp2/
- Practical Time Series Forecasting with R: A Hands-On Guid, Galit Shmueli, Kenneth C. Lichtendahl
- scikit-learn documentation: https://scikit-learn.org/stable/documentation.html
- StatsModels documentation: https://www.statsmodels.org/stable/index.html
- Prophet Python API documentation: https://facebook.github.io/prophet/docs/quick_start.html#python-api
- Forecast Package Documentation: https://pkg.robjhyndman.com/forecast/index.html
- Kaggle Dataset Hourly Energy Consumption: Over 10 years of hourly energy consumption data from PJM in Megawatts: https://www.kaggle.com/robikscube/hourly-energy-consumption
Anexo: Validación cruzada anidada para series temporales¶
En el análisis de series temporales es muy importante tener en cuenta la ordinalidad de los datos. Es precisamente ésta la propiedad característica de una serie temporal, la relación de cada punto con el todo el histórico de datos anterior a él.
Durante este proyecto, y por limitaciones de tiempo y capacidad de computación, he realizado las predicciones y los modelos con una división simple en un único conjunto de entrenamiento y un cojunto de validación. Sin embargo existe un procedimiento más sofisticado, que es la validación cruzada anidada, que podemos visualiar en la imagen siguiente:
En análisis de series temporales la proporción de datos incluída en los conjuntos de entrenamiento y validación suele ser del 80%-20%. Aunque suele haber desviaciones respecto estas cantidades dependiendo de la duración de la serie temporal, las características de estacionalidad y estacionareidad, variables covariantes, etc.
A continuación he dejado un pequeño fragmento de código que aprovecha la funcionalidad de la clase TimeSeriesSplit() del paquete scikit-learn para implementar una validación cruzada.
def plot_cv_indices(cv, X, y, ax, n_splits, lw=10, cmap_cv=plt.cm.coolwarm):
"""Create a sample plot for indices of a cross-validation object."""
# Generate the training/testing visualizations for each CV split
for ii, (tr, tt) in enumerate(cv.split(X=X)):
# Fill in indices with the training/test groups
indices = np.array([np.nan] * len(X))
indices[tt] = 1
indices[tr] = 0
# Visualize the results
ax.scatter(range(len(indices)), [ii + .5] * len(indices),
c=indices, marker='_', lw=lw, cmap=cmap_cv,
vmin=-.2, vmax=1.2)
# Formatting
yticklabels = list(range(1, n_splits+1))
ax.set(yticks=np.arange(n_splits) + .5, yticklabels=yticklabels,
xlabel='Sample index', ylabel="CV iteration",
ylim=[n_splits, -.2])
ax.set_title('{}'.format(type(cv).__name__), fontsize=15)
ax.legend([Patch(color=cmap_cv(.8)), Patch(color=cmap_cv(.02))],
['Testing set', 'Training set'], loc=(1.02, .8))
# Make the legend fit
plt.tight_layout()
return ax
# Let's separate the features_target dataframe in features matrix X and target vector y
def separate_features_target(df, target_labels):
"""
This function takes the features+target DataFrame and splits in the features and target matrices
"""
X = df.drop(target_labels, axis=1)
y = df[target_labels]
return X, y
X, y = separate_features_target(features_target, 'demand_in_MW')
display(X.head())
display(y.head())
fig, ax = plt.subplots(figsize=(7,4))
n_splits = 5 #Number of train/cv/test folds
#max_train_size = int(len(features_target)*.67/n_splits)
tscv = TimeSeriesSplit(n_splits, max_train_size=None)
for train_index, test_index in tscv.split(X):
print("TRAIN:", train_index, "TEST:", test_index)
X_train, X_test = X.iloc[train_index], X.iloc[test_index] # pandas doesn't support directly integer based slicing, so that .iloc is mandatory
y_train, y_test = y.iloc[train_index], y.iloc[test_index]
plot_cv_indices(tscv, X, y, ax, n_splits)
